str = '4phase.bmp';
global HD HR
HD = frameFilterD('linear');
HR = frameFilterR('linear');
%get input image, renormalized to [0,1]
I = readImage( str );
% I = imnoise(I,'gaussian',0,0.02);
imshow(I,[]);
%%
f = I;
psf = fspecial('gaussian',2,2);
u0 = f;
[M N] = size(f);
f = (f-min(f(:)))/(max(f(:))-min(f(:)));
usmooth =f;


%set parameter values.
global nPhases nLevels
global M N
nPhases = 4;
nLevels = 1;     %level of framelet transform
nIterations = 20;
tao = 0.4;
sig = 0.4;
gamma = 4;
beta =1;
nu = 200;
alpha = 0.5;

%  c = kmeansC(f,nPhases);
c= [0,1];
c= sort(Kmeans(I,nPhases,100));
u = cell(1,nPhases);
d = cell(1,nPhases);
for i=1:nPhases
    u{i}=f;
    d{i}= frameDec2(usmooth,nLevels,HD);
end

%set gw, the edge indicator function.
tm1 = d{1}{1};tm1(1:M,1:N) = 0;
sigma = 1;
ubar = u;
uold = u;

y = cell(1,nPhases);
for i = 1 : nPhases
    for j = 1 : nLevels
        y{i}{j}= zeros(M*numel(HD),N*numel(HD));
    end
end
wu = cell(1,nPhases);
Eold = 0;
Ecur = 0;

%BEGIN ITERATION
%-----------------------------------------------------------
for k = 1 : nIterations
    k
    tic
    for i = 1 : nPhases
        wu{i} = frameDec2(ubar{i},nLevels,HD);
        y{i} = cellPlus(y{i},cellMultiply(wu{i},sig));  
        y{i} = thresholdCell_2(y{i},1,HD);
    end



W=[];
for i = 1 : nPhases
    ry = frameRec2(y{i}, nLevels, 0, HR);
    u{i} = u{i} -tao*nu*(c(i)-f).^2 - tao*ry;
    W = [W;u{i}(:)'];
end
Wp = projSplx(W);

for i = 1 : nPhases
    u{i} = reshape(Wp(i,:),[M,N]);
end

theta = 1/sqrt(1+2*gamma*tao);
tao = theta*tao;
sig = sig/theta;
ubar = u;
%     W = [];
%     for i = 1 : nPhases
%         W = [W;u{i}(:)'];
%     end
%     [Y I ] = max(W);
%     result = zeros(M,N);
%     for i = 1 : nPhases
%         v{i} = reshape( (I==i),[M,N]);
%         result = result+v{i}*c(i);
%     end
%     figure(5);
%      colormap(gray); hold on;
%      imshow(result,[]);
uold = u;
toc
end


% % Update c
% W = [];
% for i = 1 : nPhases
%     W = [W;u{i}(:)'];
% end
% [Y I ] = max(W);
% result = zeros(M,N);
% for i = 1 : nPhases
%     v= reshape( (I==i),[M,N]);
%     count = sum(sum(v(:)>0));
%     s = sum(v(:).*f(:));
%     c(i) = s/count;
% end

W = [];
for i = 1 : nPhases
    W = [W;u{i}(:)'];
end
[Y I ] = max(W);
result = zeros(M,N);
for i = 1 : nPhases
    u{i}= reshape( (I==i),[M,N]);
    result = result+u{i}*c(i);
end
% figure;
% imshow(result,[])
%
%%
figure;
% imagesc(u0,[0, 1]); axis off; axis equal;
% colormap(gray); hold on;
% contour(u{1}, [0,0], 'r');hold on;
 imshow(result,[])









